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Abstract Hydrodynamic jets, underdense with respect to their environment by a factor of up to 10 4 , were 
computed in axisymmetry as well as in 3D. They finally reached a size of up to 220 jet radii, corresponding to 
a 100 kpc sized radio galaxy. The simulations are "bipolar", involving both jets. These are injected into a King 
type density profile with small stochastic density variations. The back-reaction of the cocoons on the beams in the 
center produces armlength asymmetries of a few percent, with the longer jets on the side with the higher average 
density. Two distinguishable bow shock phases were observed: an inner elliptical part, and a later cylindrical, 
cigar-like phase, which is known from previous simulations. The sideways motion of the inner elliptical bow shock 
part is shown to follow the law of motion for spherical blast waves also in the late phase, where the aspect ratio is 
high, with good accuracy. X-ray emission maps are calculated and the two bow shock phases are shown to appear 
as rings and elongated or elliptical regions, depending on the viewing angle. Such structures are observed in the 
X-ray data of several radio galaxies (e.g. in Abell 2052 and Hercules A), the best example being Cygnus A. In 
this case, an elliptical bow shock is infered from the observations, a jet power of 10 47 erg/s is derived, and the 
Lorentz factor can be limited to T > 10. Based on the simulation results and the comparison to the observations, 
the emission line gas producing the alignment effect in radio galaxies at high redshift is suggested to be cooled 
gas entrained over the cocoon boundary. 
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1. Introduction 

The research on simulations of extragalactic jets 
has been extended recent ly into the r e gime of very 
light j ets (compare e.g . IKrause l2002t 1 Saxton et alJ 
2002t ICarvalho fc O'Deal l2002albl : IZanni et alJ \200i 



Bicknell et al.ll2003|) . down to density contra sts (jet over 
ambient density) of rj = 10~ 5 l)Krausell2003j) . The inter- 
est has been stimulated not only by the difficulty to fix 
that parameter for real sources, but also by a first success 
in explaining some parameters of observed extragalactic 
jets. Low jet densities are needed to get large radio co- 
coon and bow shock widths. Based on the results from a 
grid of simulations and comparison to the radio and X- 
ray data of the radio gal axy Cygnus A (jCarilli fc Barthell 
Il99fil ISmith et al.ll2002lL a density contras t of roughly 
n < 10~ 3 has been claimed for this source ijRosen et alJ 
Il999t lKraus"ell2003l) . All these simulations have been car- 
ried out in axisymmetry and a scale of a few dozen jet 
radii. The aim of this paper is to verify and extend the 



earlier simulation results on the 100 kpc (=200 jet radii, 
R\ ) scale a n d in t hree dimensions. 

IKrause! l|2003() found that very light jets start their 
life in a blastwave-like phase. During that phase, termi- 
nating when the bow shock reaches R\ ~ Rj/2rj 1 / 4 , the 
bow shock is spherical because of pressure predominance, 
obeying the blastwave equation of motion: 



M(r')r' dr' = 2 / dt' / E(t")dt" 



(1) 



Here, A4(r) is an arbitrary spherically symmetric ambi- 
ent gas mass distribution and E(t) is the energy injection 
law. For stationary energy injection E = Lt and constant 
ambient density po, one gets: 



5Lt 3 



1/5 



(2) 
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This equation relates the total luminosity, L, of the jet 
to observables like the bow shock radius and its veloc- 
ity (via time derivation of (0). With information on the 
bow shock propagation at later phases, it can be hoped 
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that the jet power could be reconstructed from the bow 
shock shape, which is probably already observed in the 
case of Cygnus A by Chandra (see sect. |5J. It is there- 
fore important to check the evolution of the bow shock 
at later phases. This was difficult for most simulations 
so far because of the assumed reflection symmetry in the 
equatorial plane of the system. Interaction of the back- 
flow with this boundary disturbs the sideways evolution 
of the bow shock. Therefore, the symmetry assumption 
has been dropped in the present work, and the simula- 
tions are bipolar, evolving back-to-back jets in opposite 
directions, a technique that has emerged only recently 
llRevnolds et all 120011 1 Krause fc Camenzindll2002l l2003t 
iBasson fc Alexanderll2003|) . A fully 3D simulation is pre- 
sented in sect. |3 and a large scale axisymmetric simula- 
tion is presented in sect. 0] Comparisons to observations 
are presented in sect. [5] Part of this is an update of pre- 
viously published results, using a more involved analysis. 



2. Numerics 

The magn eto-hydrodynamic (MH D) code NIRVANA was 
employed i Ziegler fc Yorkelll997l) . It solves the hydrody- 
namic equations in three dimensions (3D) for density p, 
velocity v, and internal energy e: 



"ext ext "ext 



dpv _ . . 
+ V • (pvv) 



= 



-Vp - 

-p V • v - p 2 K 



(3) 
(4) 
(5) 



where $ denotes an external gravitational potential and 
A is a cooling function. Here, bremsstrahlung (A = 7.5 x 
10 20 VT(1. + 4.4 x KT 10 T) cm 3 crg/s) has been used for 
the 3D computation. 

The code was vectorised and parallelised by OpcnMP 
lik e methods, and successfully tested on the NEC SX- 
5 ijKrause fc Camenzindl 120021) at the high performance 
computing center in Stuttgart (Germany), where the com- 
putations have been carried out. 



2.1. Boundary conditions for the 3D cylindrical grid 

For the 3D simulation, a cylindrical grid was employed. 
The disadvantage of the cylindrical coordinates (Z, R, <j>) 
compared to the Cartesian ones is the appearance of 
internal boundaries. The grid somehow has to be con- 
nected across these boundaries. In the <fi direction, peri- 
odic boundary conditions were applied. For the boundary 
on the axis, no analogue could be found in the literature. 
The boundary condition here is similar to the periodic 
case: One side of the grid should know about the other 
side. Therefore, three cells were used below the axis, to 
which consequently the index ir ■ = 3 was assigned. With 
this choice and the staggered mesh, equations I ilM are 
well defined everywhere on the grid. For the scalar quan- 




Figure 1. Sketch of the cylindrical grid. 



tities, the following boundary conditions were applied: 

f(iz,iR,U) ■= f(iz,6 — iR,i<i> + iir) ,ir = 0,1,2 
f(iz,3,U) := < f{iz,Z,i<i>) > | lz =const , 

where denotes the number of grid cells corresponding 
to the angle tt. The second equation indicates that di- 
rectly on the axis , the scalars were averaged over for 
constant iz, since all the different refer to the same 
physical point. Due to the staggered mesh, the Z and <fi 
components of vectors are shifted by half a grid cell away 
from the axis. Therefore for these components, there arise 
no special problems, and the boundary condition is: 

Vz,<j>(iz,iR, i<j>) ■= vz, <p(iz, 5 - i R , i^ + i^) , ir = 0, 1, 2 

The radial vector components are not shifted away from 
the axis. So in principal, they all denote the same physical 
point. However, the flow should be allowed to cross the 
axis from one side to the other. This is only possible if the 
radial velocity takes a reasonable value there. Therefore, 
two possibilities arise for the boundary conditions: 



VR{iz,iR,i<i>) 
v R (iz, 3, 14,) 



VR(iz,6~-iR,i<f> + in) ,tR = 0,1,2 

: I 

{vr{izAiH) ~ Vr(iz,2, v))/2 : II 



In order to check the influence of these two different 
boundary conditions on the simulation, the 3D jet was 
simulated twice, the first time with the case I, and the 
second time with the case II boundary condition. For both 
runs, the integrated quantities (directionally split energy 
and momentum, and mass) and also the timesteps were 
equal. Also, from comparison of the contour plots, no dif- 
ference could be found. Hence, the flow seems to have 
enough possibility to flow past the axis, and the detailed 
behaviour at that line does not influence the result by 
much. The simulation result also shows that this approach 
in general works in regions of undisturbed jet flow. Where 
the jet is dominated by instabilities and wants to bend, the 
axis looks like an obstacle. Therefore, the method seems to 
be problematic if details about the jet beam are of special 
interest. Nevertheless, the results confirm that the rep- 
resentation of the beam is acceptable, and for the other 
regions the output was fine. The big advantage of the 
cylindrical grid is the reduction of the number of cells 
to compute. The reduction of the physical computational 
volume is 22%. Since the ^-direction can be covered with 
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fewer cells than a third dimension in Cartesian coordi- 
nates, this approach saves a significant amount of mem- 
ory and CPU time. The computation would not have been 
possible on a 3D Cartesian grid of similar central resolu- 
tion. 

3. 3D simulation 

3.1. Simulation setup 

A cylindrical grid was used for the jet simulation (compare 
section EU). The size of the computational domain was: 
Z G [-69 kpc,69 kpc], R € [0,57 kpc] and <f> G [0,2tt]. 
2042, 805, and 57 grid points were used in the Z,R and 
4> directions, respectively. With a jet radius of r- s = 0.55 
kpc, this gives a resolution of 8 points per beam radius 
(ppb). The resolution in 4> direction scales linearly down 
from 16 ppb on the jet boundary to 0.2 ppb on the edge of 
the grid. The grid was initialised with an isothermal King 
cluster atmosphere: 

p c (r)= Pc .oi [ l + ^j , (6) 

where r = \/ R 2 + Z 2 denotes the spherical radius, p c .o = 
1.2 x 10~ 25 g/cm 3 is the characteristic density, (3 = 0.75 
and a = 35 kpc is the core radius. In order to break the 
bipolar and axial symmetry, this density profile was mod- 
ified by random perturbations of the following kind: 

1. With 10% probability the density was increased by a 
factor between 1 and 1.4. 

2. With 5% probability the density was increased by a 
factor between 1 and 2, only if the cell was unmodified 
in the first process and the Z coordinate was positive. 
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Figure 2. Bow shock shapes for the 3D run at t = 2.04 Myr. 
On the top panel the shape in one meridional plane is compared 
to the elliptical R = 7.5^/1 - (z + 1) 2 /169. The bottom panel 
compares the bow shock shape for different meridional planes 
to the parabola R = 2.3(26- Z) 1/3 (thick line). The bow shock 
is axisymmetric in the middle, where it can be well represented 
by an ellipse. For \Z\ > 10, the bow shock is bumpy and not 
axisymmetric. That part is called cigar-like. 
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Figure 3. Radius where the bow shock hits the line (</> — 
0, Z = 0) versus time, with a globally fitting function. 



The jet is injected in the middle of the grid in the re- 
gion Z G [-0.55,0.55 kpc], R G [0,0.55 kpc], and 
4> G [0,27r]. This region has the constant values: pj c t = 
6.68 x 10 - 28 g/cm 3 , Vz = ±0.4c, where c denotes the 
speed of light. The plus sign applies for the positive Z re- 
gion and the minus sign for the negative one. The kinetic 
jet luminosity is Lkin = 1-04 x 10 46 erg/s for both jets to- 
gether. The pressure was set in order to match the exter- 
nal pressure at that position. This gives a slightly varying 
density contrast across the grid of r\ = pj e t/ Pcxt ~ 7x 10 -3 
and an internal Mach number M = 10. The jacket of the 
injection cylinder is a further boundary. This boundary 
is open, so the material can leave the grid here into the 
central kiloparsec of the simulated radio galaxy, which is 
not attempted to model here. The temperature in the ex- 
ternal medium is set to T = 3 x 10 7 K. The cooling time 
in the shocked cluster gas is approximately 100 Myr. The 
jet is expected to propagate through the whole volume in 
10 Myr. So, cooling by bremsstrahlung marginally influ- 
ences the state of the gas. This was taken into account 
(comp. ||SJ and sect.0). Thus the calculation is not scal- 
able anymore, formally. But given the smallness of the 
effect, scaling should be possible, in practice. Because the 
atmosphere is isothermal, the pressure varies in the same 
way as the density. In order to keep the system in hydro- 
static equilibrium, gravity by an assumed dark matter dis- 
tribution had to be applied. The gravitational potential, 
necessary to prevent the King atmosphere from exploding 
is: 

$ DM = M^log(l + (r/a) 2 ) , (7 ) 
2/imn 

where p is the number of particles per proton mass. 

3.2. Early evolution &. verification of boundary 
conditions 

A testrun was performed on a ten times smaller grid, in 
order to verify the boundary conditions. This early evolu- 
tion is shown in Fig. The images show the formation of 
a bow shock and backflows, for both jets. 
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Figure 4. The very early evolution of the 3D simulation. The meridional slices at <f> — 0, tv, joined at the axis, of number density 
(left) and Mach number (right) are shown at 0.03 Myr (top) and 0.06 Myr (bottom). Values below 0.1 were omitted in the plot 
of the Mach number. The bottom row shows that sometimes small numerical artefacts at the axis are present, especially in the 
Mach number, which includes the radial velocity. 
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Figure 5. Slices of the tracer at t = 0.32 Myr for different axial positions. The shocked ambient gas is shown dark. Compressed 
regions are darker, rarefied and diluted regions appear lighter. Cells that contain less then roughly 10% ambient gas, i.e the 
beam plasma, are shown in white. 



In that phase, the jet almost ignores the stochastic na- 
ture of the density. The bow shock has a round and regular 
shape, and the density varies smoothly along its surface. 
However, the evolution on the two sides is different. The 
jet on the left-hand side, where the average background 



density is smaller, evolves faster: It produces a bigger bow 
shock, and a faster backflow at equal times. At t = 0.06 
Myr, the two bow shocks are nearly joined together. As 
the later evolution shows, a single round bubble forms, 
soon after. The upper and lower halfs of the pictures fit 
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Figure 6. Meridional slices for <j> = 0/n at 0.32 Myr. The quantity plotted is indicated on top of the individual figures, units 
can be found next to the bar. The radial Mach number is defined to be positive away from the Z-axis. A positive toroidal Mach 
number is intended to mean motion into the plane for r > 0, and out of the plane for r < 0. The beam shows rotation in some 
places and translation in some others. Jet beam material can be decerned by its high entropy index. 



quite good together, although there is sometimes a suspect 
spike at the apex of the bow shock on both sides. This nu- 
merical artefact reminds us about the imperfect treatment 
of the axis. The effect is rather small, and even absent for 
most of the simulation time. This supports the choice of 
boundary conditions on the axis (compare section fe.ljl . 

3.3. Medium term evolution 

A snapshot of several quantities at t = 0.32 Myr is shown 
in Fig. |5| The jet plasma takes the lowest density, and 
the highest temperature values. At that time, the evolu- 
tion continued in the same way as in the early phase: the 



right-hand side, propagating into the on average denser 
medium, develops a broader cocoon, and is slower. The 
once separated bow shocks have united. The shape of the 
bow shock in that phase is oval, a sign of the blastwave 
phase. The bow shock shows two extensions in jet direc- 
tion. Due to their later appearance, these parts of the 
bow shock will be called cigar-like. The aspect ratio of 
the bow shock is nearly 2, thanks to the extensions, which 
contribute approximately 0.5 to the aspect ratio at that 
time. 

The plots of the Mach number show a slightly super- 
sonic backflow. The radial Mach number is positive for 
motion away from the axis, the toroidal Mach number 
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Figure 7. The same as Fig. [(J but for t = 2.04. (Only some of the variables are shown. 



changes sign at R = for motion out of (into) the plane. 
The flow in R and Z directions is well ordered, whereas 
in the azimuthal direction, arbitrary fluctuations are seen. 
The jet beam shows both rotation (same colour above and 
below R — 0) and translation away from the axis (colour 
change at R = 0) . The latter supports the chosen approach 
as being able to represent beam motions away from the 
axis of symmetry. 

The pressure plot shows a high pressure at shocks in 
the beam, especially on the axis. While this is in principal 
correct, the exact amount of the pressure could be influ- 



enced by the choice of the cylindrical coordinate system. 
The first shock in the beam is stronger on the left-hand 
side. This follows from the higher pressure there, but also 
from the higher inclination angle to the axis. Since the 
beams have identical conditions, why are the shocks not 
identical? The shocks in the beam are driven by what- 
ever hits it. From the density plot, this seems to be the 
cocoon on the left-hand side, and the entrained cluster 
gas on the right-hand side. The backflows collide approx- 
imately at the center, forming a region of enhanced pres- 
sure. This region is asymmetric at the time shown here. 
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Due to the stronger backflow from the right-hand side, the 
region has moved to the left. The higher pressure there 
drives a stronger shock into the beam. This explains why 
the first shocks arc different on both sides. 

A passive tracer variable was advected with the flow, 
set initially to unity in the shocked ambient gas, to zero 
in the right beam, and to —1 in the left one. This tracer 
enables differentiation between the beam and the ambi- 
ent gas. Slices at constant axial position of the tracer 
at t = 0.32 Myr are shown in Fig. The right side 
has an approximately axisymmetric cocoon, while the left 
one is clearly not (compare slices at Z = —3.5 kpc; 
Z = —1.5 kpc). The slice at Z = —0.83 kpc demonstrates 
how the cocoons manage to merge together: The left one 
is smaller and slips inside the right one. This also brings 
shocked ambient gas inside the cocoon, which can still be 
spotted in the Z = 0.5 kpc slice. 

The stronger shock on the axis on the left-hand side 
causes the beam to widen, due to high pressure. On the 
right-hand side, less energy has been converted into heat, 
the beam stays narrow and delivers more power to the 
terminal shock, where the pressure is consequently higher. 

Fingers of shocked cluster gas, reaching down to the 
beam surface, are present just as t hey are in the unipo- 
lar 2 .5D simulations (compare e.g. iKrause fc Camenzindl 
l200l|) . They are generated in the following way: Kelvin- 
Hclmholtz instabilities are excited at the boundary be- 
tween cocoon and shocked external medium. The back- 
flow advects those instabilities, while they are amplified. 
Hence, they are biggest in the symmetry plane. The mate- 
rial in the symmetry plane could in principal flow outward, 
creating bumps in the bow shock, or inward. The simu- 
lation clearly shows no sign of outward motion. Instead, 
the gas is transported far down into the radio cocoon in 
geometrically thin fingers. Soon, their extension falls be- 
low the resolution limit of the simulation, and the gas 
mixes with cocoon plasma. However, in reality, the two 
phases may remain separate. The magnetic field of the 
radio plasma further supports the separation of the two 
phases. In principle, in that way new fuel could be chan- 
neled to the central source. However, the dynamics of the 
central kpc is beyond the scope of this work. The central 
high pressure region pushes the entrained shocked cluster 
gas to the right, where it slips between beam and cocoon, 
thereby widening the right cocoon. 

It has already been pointed out above that at the time 
the simulation is shown in Fig. H3 0.32 Myr, the left jet 
converts more kinetic energy into heat than the jet on the 
right-hand side. At that time the tip of the jet has al- 
ready advanced approximately 20 % further towards the 
left than towards the right. The difference of the average 
density is only 2.5 %, which cannot explain the fast prop- 
agation of the left jet, considering the estimate by the one 
dimensional force balance: Vhead ~ -y/^evbeam, where e is 
the ratio of beam to head area. One important result is 
therefore that the nonlinear dynamics amplifies the effect 
of the different density on both sides by more than a factor 
often. The situation at that time is unstable, and the later 



timesteps show that the right jet catches up, and outruns 
the left one not later than at t = 0.95 Myr. The situation 
stays that way until the end of the simulation. On average, 
the right jet is approximately 10% faster than the left one, 
at late times. This is in conflict with naive intuition, but is 
readily explained, considering that the stronger backflow 
from the right jet shifts the central pressure enhancement 
to the left, where stronger oblique shocks are driven into 
the beam slowing the left jet down, as pointed out above. 

3.4. Long term evolution 

The final snapshot at t = 2.04 Myr is shown in Fig.0 The 
plots show the usual picture of a hydrodynamic jet sim- 
ulation, largely consistent with FR II radio galaxies: The 
cocoon is now nicely placed around the jet beam. The as- 
pect ratio of the bow shock is 3.6. The Kelvin-Hclmholtz 
instabilities show up prominently. They still grow towards 
the center, and develop into long fingers at the innermost 
positions. The pressure shows a regular spacing of shock 
compression and rarefaction zones in the beam. High pres- 
sure regions are small and show up only at the end of 
the beams, where the Mach disk is located. The oblique 
shocks in the beam are now weaker than at t = 0.32 Myr. 
This follows from the smaller angle with the jet axis. The 
central region, with a diameter of roughly 10 kpc, is now 
dynamically calm. No large Mach numbers are observed 
there, and the pressure is approximately constant. The 
density takes intermediate values. This is now a relaxed 
region where jet plasma and shocked cluster gas are mixed 
(mixing may not happen in nature, see above). 

3.4.1. The shape of the bow shock 

Figure [2 shows the bow shock shape for the final snapshot 
in detail. It has an axisymmetric part in the middle, where 
it can be well represented by an ellipse. The center of 
this ellipse is not the origin of the grid, but is shifted 
to the left by one kpc. Such a shift of the center is also 
found in the larger 2.5D simulation. The elliptical shape 
ends at \Z\ f=a 10 kpc where two cigar-like extensions join 
the bow shock. These extensions are 3D in nature with 
several bumps. They can be represented, on average, by a 
parabola of rank three. 

3.4.2. The law of motion of the bow shock 

In all the plots in the long term evolution, the two bow 
shock phases are easily discernable. The inner bow shock 
structure is a remnant from the blastwave phase. In its 
early phase, the radius of the bow shock should have 
obeyed equation (J2J), i.e. r oc t 6 . (This part of the bow 
shock is far inside the core radius wherefore the density 
profile is roughly constant.) In the following, this law is 
checked for early and later evolution using the bow shock 
position in the Z = plane. For every hundredth timestep, 
the positive radial bow shock position for <f> = and Z = 
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was saved, up to the final timestep 220200, corresponding 
to t = 2.04 Myr. These values are plotted against time 
in Fig. [3J The positions were fitted with a function of the 
type: a + bt c . The constant a was allowed for in order 
to take into account an artificial offset due to the ini- 
tial conditions imposed by a jet with a final radius and 
non-developed structure at t = 0. The fit was carried out 
globally, and separately for the early and the late evolu- 
tion. The result is given in Table ^ The global fit gives 

Table 1. Fit parameters for the bow shock position. The 
star denotes a fit with fixed a = 0. 



time range [Myr] 


a 


b 


c 


global 


1.02 


4.44 


0.58 


[0 : 0.5] 


0.41 


4.52 


0.38 


[1:2] 


1.25 


4.21 


0.60 


[1 : 2]* 





5.45 


0.49 



3.4.3. Implications on the jet power 

From the b parameter of the late evolution, one can cal- 
culate the jet power, according to equation lf2"|: 

L kin = 2.23 x 10 67 po& 5 w 3.75 x 10 45 erg. (10) 

This is approximately 40% of the true jet power (it would 
be 50% taking b from the global fit, and 130% for the late 
fit with a = 0). One can also use the bow shock position 
and velocity at e.g. t = 2.04 Myr in order to compute its 
power. From one gets: 

L = 1007rp (r bow ) 2 (w bow /3) 3 . (11) 

This yields a jet power of L m 10 46 erg/s, accurate to 
a factor of a few because of the slow bow shock veloc- 
ity and the low spatial resolution. The true jet power is 
1.04 x 10 46 erg/s. It follows that the information on the 
jet power is conserved in the law of motion of the central 
part of the bow shock at least with an accuracy of a factor 
of a few. 



an exponent c of 0.58, quite close to the exponent analyt- 
ically derived for the blastwave phase (0.6). However, the 
behaviour is quite different at late times compared to the 
early evolution. At late times, the exponent is in the range 
0.5 to 0.6 (compare below), whereas for the early evolu- 
tion, c is close to 0.4, which would be the value expected 
for a supernova blastwave (initial energy input with no 
additional supply, compare Q): 



15£ t 2 
4-irpo 



1/5 



(8) 



This is due to the initial conditions of the simulation: in 
order to get a propagating jet, one has to start with a 
short propagating jet. The energy of that initial jet is quite 
high. It can be estimated by calculating the energy stored 
in the initial injection box: Eq ps irR?hpjV? , where h is 
the height of the the cylindrical region. This amounts to 
approximately 10 57 erg. The energy is also given by the 
parameter 6, according to equation ©, and since the den- 
sity is known: 



E = 2.35 x 10 8U po^ ~ W erg. 



(9) 



The difference of a factor of ten could be attributed to in- 
creased heat production in the very early evolution, where 
the unphysical initial condition, which is always present in 
such simulations, relaxes into a regular jet structure. 

At late times, the initial conditions may be regarded 
as relaxed. In that case, the fit parameter a could be set 
to zero. This would change the best fitting exponent from 
0.6 to 0.5. An exponent lower than 0.6 would be expected 
due to the increase of the aspect ratio. However, the effect 
is quite small. 



3.4.4. Emission maps 

The emiss ion due to optically thin brcmsstrahlung (e.g. 

was computed, mapped onto a Cartesian grid, 
and integrated for different viewing angles (Figs.|Hland|njl- 
The general X-ray emission properties of shocked 
ambient gas in jet simulations have been discussed by 
Clarke et a which has been updated recently by 



Zanni et alJ l|2003l) . The idea is that the gas is pushed 



aside by the jet cocoon. Depending on its compression, it 
will or will not form X-ray deficits at the location of the 
cocoon, and bright shells at the edges. Averaging the gas 
distribution and neglecting flows in the shocked ambient 
gas, the critical parameter is the relative shell thickness 
£, defined as the width of the shocked ambient gas region 
divided by the local bow shock radius. Then, the ratio of 
observed flux to the flux of the initial condition is given 
by: 



/7/=^r 1 (2-0" 



(12) 



where T depends on the pre-shock and post-shock temper- 
atures, and is close to unity. Therefore, X-ray deficits could 
be observed for a shell thickness £ above 38%, for sources 
located at 90° inclination. Here, the shell thickness is com- 
paratively low for the whole simulation. Therefore, at high 
inclinations the X-ray surface brightness never falls below 
that of the undisturbed King atmosphere. This picture 
implies that the deficit is pronounced for low inclinations, 
since most of the gas is shifted aside. Indeed, wc find these 
deficit regions in Fig. for angles of 10° and lower. The 
deficit is much more evident at the later time. This is due 
to the different morphology: At t = 0.32 Myr (Fig. El the 
bow shock shape is essentially elliptical with only a small 
cigar-like extension. In that phase, the gas is compressed 
more isotropically (compare Fig. 0) . The thickness of the 
shocked ambient gas shell is roughly the same in all direc- 
tions. At the later time, the prominent cigar has displaced 
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Figure 8. Bremsstrahlung emission maps for the 3D run at t = 0.32 Myr. From top to bottom, the viewing angle is 0°, 10°, 
30°, 60°, and 90°. The left column shows the emission map, the middle one a vertical slice through the center, and the right 
one a horizontal slice through the center of the emission map. For comparison, the undisturbed cluster emission is given. The 
bow shock morphology is reflected in the ring like structures of differing sizes and surface brightness profiles. 



the gas sideways, producing the deficit when viewed pole 
on. Another reason is the decreasing density. The region 
in front of the jet head is furthermost from the center and 
the gas is therefore thinnest. 

The two phases of the bow shock, cigar and elliptical 
(comp. sect. I3.4~TI . show up prominently in the emission 



maps. They form circular and elliptical rings, respectively, 
depending on the viewing angle. Where the rings partially 
overlap, they are brighter, producing the impression of 
ring segments (e.g. Fig. El 10° and 30°., Fig.^J 10°). The 
structures can also intersect on the line of sight, producing 
bright spots (Fig. EI 30°). The pole-on figures show at least 
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two rings: one from the cigar phase and one from the inner 
elliptical bow shock part. 



Initially, the edges arc brightened for any viewing an- 
gle. This changes at late times for high inclination. Here, 
the edges are prominently brightened only for the direc- 
tion perpendicular to the jet axis. The reason is the declin- 
ing density profile that is superimposed on the jet features. 
Edges and rings are typically enhanced by a factor of a few 
in surface brightness. They should therefore be detectable 
for some sources in this phase. 



4. 2.5D simulation 

4.1. Simulation Setup 

In order to study the jet evolution on larger scale, a 2.5D 
simulation was performed with similar initial conditions to 
the 3D simulation in the previous section. The simulation 
was run for 20 Myrs of simulation time; during that time 
the jet reached an extent of 110 kpc which corresponds 
to 220 jet radii. In the radial direction, the jet reached 
31 kpc. The jet radius was set to 0.5 kpc and the resolu- 
tion was set to 20 cells per jet radius. This corresponds 
to [4400 x 1240] cells in total. With that resolution global 
parameters like the bow shock velocity on the Z-axis or en- 
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ergy and momentum conserv ation are accurate to « 10% 
ijKrause fc Camenzindll200l|) . p c $, a and (3 were chosen 
slightly different to the previous model, because the side- 
ways expansion of the bow shock should be able to escape 
the constant density part during the simulation time. The 



Bow shock shape at t= 1 Myr compared to fit 



values arc: p c ,o = m p /cm 



lOkpc and /? = 0.75. The 



temperature was again set to T = 3 x 10 7 K. The jet is 
injected with a density of pj c t = 10~ 4 x p c , the sound 
speed in the jet was set to 20%c, and the jet's Mach 
number to M = 3. The high internal sound speed and the 
low Mach num ber was chosen because a parameter study 
llKrausell2003T) has shown that very light jets adjust their 
pressure quickly to the cocoon pressure by oblique shocks 
(no expansion or contraction) and that high Mach num- 
bers cannot be sustained in a non-relativistic very light jet. 
No cooling was taken into acc ount, and hence th e simu- 
lation is scalable as outlined in IZanni' et all fcOO.I) . With 
these parameters, the simulation took about a month on 
an SX-5 supercomputer. It was therefore not possible to 
do this in 3D, for now. 



4.2. Results 

We present logarithmic density plots of the simulation re- 
sults for four different simulation times (5, 10, 15, 20 Myr) 
in Fig. 1111 The morphology that appears in these figures 
is a continuation of previous simulations that could not 
reach the size shown her e. The 5 Myr figure shows the 
state that was reached bv lKrausel 1 2003h . and extensively 
discussed therein. In this early phase, the bow shock is 
spherical, its radius following an expansion law given by 
the force balance equation 0J. At later times the cocoon 
transforms via a conical phase towards a cylindrical one. 
This is a remarkable result because classical double radio 
galaxies also often have cylindrical cocoons. 



4.2.1. The shape of the bow shock 

The bow shock surface was extracted and fitted by ellip- 
tical functions as in the previous section. The result is 
shown in Fig^| As in the 3D case, the bow shock can be 
represented by an ellipse in the center and by a parabola 
near the tip. This parabola now has rank two while it was 
a rank three parabola in the 3D case. The difference is 
due to the different jet nature. The beam in the 3D simu- 
lation stays intact all the way to the tip of the bow shock. 
Therefore it can deliver its momentum to a smaller area. 
In the 2.5D simulation the beam becomes turbulent, en- 
training cocoon and shocked ambient gas, delivering the 
momentum to a larger area. Also here, the center of the 
ellipse is shifted to the right, by 1.3 jet radii. This corre- 
sponds to the armlength asymmetry, i.e. the right jet is 
longer for t > 3 Myr, by a few percent. This confirms the 
findings of the 3D simulation. 




Bow shoe* shape at t= 20 Myr cornparsd to fits 




Figure 10. Shape of the bow shock at 1 Myr (top) and at 
20 Myr (bottom) for the 2.5D simulation. The fit functions 
are: R = 6.61 x ^1 - (Z + 0.23) 2 /53.28 (1 Myr), and R = 
3. 66 x VZ + 55.6, R = 3.91 x V55.5 - Z, and R = 30.7 x 
y/l - (Z + 0.74) 2 ./1277 (20 Myr). This means that the bow 
shock evolves from almost elliptical to elliptical + parabola 
extensions. 



4.2.2. Beam stability 

The beam in the 3D simulation was stable, whereas the 
beam in the 2.5D simulation becomes turbulent. The criti- 
cal factor for this behaviour is not so much the dimension- 
ality (that should work in the opposite direction, although 
the stability is probably somewhat enhanced by the choice 
of the cylindrical coordinates) but the density contrast. 
Lighter jets a utomatically moderate their Mach number 
llKrausdbool . This can be understood from the pressure 
equilibrium within the whole jet. Because in a very light 
jet, the sound speed is high and the expansion speed is low, 
the pressure is almost equal anywhere in the jet (compare 
Fig. 112(1 . The jet beam adjusts to that pressure by varying 
the strength of its oblique shocks. At the actual location 
of the shocks, the pressure is higher. This sums up to a 
decline according to a power law with exponent —8 in the 
histogram and can be seen in the pressure distribution 
(Fig. 112(1 . The sound speed in the beam is therefore above 
10% the speed of light, which makes it hard for the jet to 
be highly supersonic in a non-relativistic simulation (the 
typical Mach number in the beam is between one and two) . 
The problem is even more severe at earlier times because 
the average pressure decreases with time (see below) . It is 
well-known that jet beams at low Mach number are dis- 
rupted quickly (e.g. iBodo et al.lll994ll995tll99s|) . 

4.2.3. Pressure evolution 

The average jet pressure is the driving force of the in- 
ner elliptically shaped part of the bow shock. Figure IT2"1 
shows that the pressure in the jet system monotonically 
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Figure 11. Four snapshots of the 2.5D simulation. The logarithm of the number density is shown. The times for the snapshots 
are indicated on top of the individual figures. The same part of the grid is shown in each case. The jet forms first a spherical 
bow shock, associated with a spherical cocoon. The cocoon then transforms via a conical state to a cylindrical one. The bow 
shock's aspect ratio (length/width) also grows, up to 1.8. 
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Figure 12. Left: Pressure distribution for the 2.5D simulation after 20 Myr. Top right: Axial pressure average over the radius. 
Bottom right: Pressure histogram within the jet (i.e. beam, cocoon, and shocked ambient gas). It can be represented by a broken 
power law with humps around the most frequent value. The power law indices are 5 and —8, respectively. 90% of the volume 
has a pressure in the range [1 — 2] x 10~ 9 dyn/cm 2 . The noise represents the statistical error. 



decreases with radius. Close to the axis, the pressure 
is higher because of shocks in the beam region. In the 
shocked ambient gas region, a new equilibrium of grav- 
ity and pressure appears. The smallest pressure values are 
located at the bow shock, roughly 20% below the average. 

In the previous section, it was shown that the inner 
part of the bow shock roughly follows the spherical prop- 
agation law inside the core radius. The accuracy of 
this law will be checked in the following also for the larger 
2.5D simulation. In this case, the bow shock has propa- 
gated more than three core radii in the sideways direc- 
tion. Notice that the formulae derived for blastwaves with 
strong shocks arc still applicable here, as shown in the ap- 



pendix. In the spherical approximation, the average pres- 
sure inside the whole jet is given by : 

Pj = (7-1) y , (13) 

Here, 7 = 5/3 is the adiabatic index, and Vj is the jet vol- 
ume (including beam, cocoon and shocked ambient gas). 
The power L includes all sources of energy, i.e. the flux 
of kinetic and internal energy through the jet channel, 
the flux of internal energy entering through the surface 
of the bow shock, and the energy lost by work against 

1 The formula neglects the increasing but small part of ki- 
netic energy stored in the beam. 
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Figure 13. Approximate internal energy flux through the bow 
shock (Lb,mt) and energy loss rate due to gravitational uplift- 
ing (Lgrav) over simulation time. Lmt.c = 4-7rr p e u/(7 — 1) was 
calculated using the bow shock's radius and sideways veloc- 
ity from the simulation. The fitted functions are: 0.57 1 37 (0- 
5 Myr), and 1.51 i~°' 08 (10-20 Myr). The gravitational energy 
loss rate levels off at 30% of the beam energy flux. 
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Figure 14. Average jet pressure over simulation time. The 
stars show the values measured in the simulation, crosses 
mark the pressure according to a spherical approximation, 
plusses show the relative difference between the former (in 
%). Corresponding symbols are connected with solid lines. The 
lines show fits to the pressure measured from the simulation. 
The fits are: 32.84i _069 (0-5 Myr, solid line), and 227.95 t~ lm 
(15-20 Myr, dashed line). The best fit for the spherical approx- 
imation in the range 15-20 Myr is: 126.05 1 -1,50 (not shown). 



the gravitational field. The simulation takes these effects 
properly into account. In order to test the pressure evolu- 
tion against deviations from the simple spherical law these 
power sources are calculated as follows: 

Via the jet beam, a constant power is injected, given 
by: Lj = Lj.kin + ij.int = nRfpvf(l + (7(7- l)M 2 /2)~ l ) = 
9 x 10 45 crg/s. 



Through the jet's bow shock, the power Lb.int = 
47rr^WbPKing(?'b)/(7— 1) enters the jet system in the spher- 
ical approximation, and are the bow shock's po- 
sition and velocity, which is measured at the position 
of the greatest cylindrical radius. PKing denotes the ini- 
tial pressure distribution of the isothermal King profile. 
In the spherical approximation, Lb.int should be propor- 
tional to t s ( K +3)-i^ w here k would be the exponent of an 
isothermal power law distribution in pressure and den- 
sity (p, p cx r K ), and 5 the exponent of a local power 
law approximation to the sideways bow shock propaga- 
tion (r oc t 5 , 5 = 3/(k + 5)). But for a King profile, it 
takes some time until 8 adjusts to the value it should have 
due to the local power law approximation. The latency 
is caused by the memory of the flow of its history due 
to the fixed amount of energy and mass at a given time. 
£b,int is shown in Fig 1131 together with fits at early and 
late times. It should first rise proportional to t 4 / 5 . The fit 
gives an exponent of 0.37 which again reflects the effect 
of the initial conditions. The fitted exponent for Lb,int of 
—0.08 after 15 Myr is in agreement with a local power law 
exponent for the pressure of k = —1.84, where 5 = 0.8 has 
been adopted from the measured bow shock propagation 
for that time span. The local k in the fitted region is in 
the range -1.7 to -2.0. 

The rate of change of gravitational energy, L gra v, is 
computed from the simulation data, averaging over one 
million years and is plotted in Fig^J L gISV rises from one 
to thirty percent of the jet beam's energy flux, where it 
seems to converge towards the end of the simulation. The 
gravitational energy loss rate and the power entering as 
internal energy through the bow shock make up a similar 
contribution to the energy within the jet system as the 
energy flux through the beam. 

Using L = Lj + Lb,int - L gISV , where £ b ,int and L grav 
now denote the time-averaged value at a given time, H13JI 
can be evaluated, where v and t are given by (Q: 



Lt 2 
Mr 



- / M(r')dr' 
L Jo 



1/3 



(14) 
(15) 



and r denotes the sideways extent of the bow shock. 
FigurelT4l shows this analytical estimate together with the 
data from the simulation. Here, r was related to time via 
measurement from the simulation. The agreement is quite 
good, in general. The analytical formula follows the slope 
of the simulation data, but underestimates it by up to 
« 20%. 

Towards the end of the simulation, the pressure evolu- 
tion can be approximated by a power law with exponent 

— 1.66. For the same time range, the spherical approxima- 
tion is well approximated by a power law with exponent 

— 1.50. Note that these are only local approximations to 
curved functions that have not yet reached the asymptotic 
power law regime. 
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Figure 15. Bow shock radius at Z—0 versus time (squares) 
with fits and compared to spherical approximation, includ- 
ing all power sources (see text, plus-signs). The three fits 
are: 4.66 i 53 (0-5 Myr), 2.99 1 ' 78 (15-20 Myr), 2.36 1 ' 81 (15- 
20 Myr, spherical approximation). 

4.2.4. Sideways motion of the inner bow shock part 

From the pressure evolution (previous section) one should 
expect that the sideways motion of the inner bow shock 
part roughly follows the spherical approximation for the 
whole simulation time. The density profile used here has 
the asymptotic power law approximations: lim ri _>o(p) = 
Po, and lim rt _ +oc (p) oc r~ 9 / 4 . Therefore, the bow shock 
should expand with r oc £ 6 at the beginning, steepening 
towards r oc £ 109 , at least as long as it remains spheri- 
cal. The radial bow shock position was determined every 
0.4 Myr (Fig. 115(1 . The resulting function has a curved 
nature. Following the arguments of sect. 13.4.21 the bow 
shock propagation was locally fitted by a function of type 
a + bt c . The resulting parameters for the different regions 
are shown in Tabled Usually, a — 0, since only the very 
late evolution of the jet is studied. Only for the time span 
up to 5 Myr a fit with a =/= has been included because 
here it is possible that effects from the initial condition 
still dominate the propagation. For comparison, also fits 

Table 2. Fit parameters for the bow shock position in the 
2.5D simulation. The star denotes a fit with fixed a = 0. 
Two stars denote fits to the spherical approximation with 
fixed a = 0. 



time range [Myr] 


a 


b 


c 


[0:5] 


2.00 


2.67 


0.76 


[0 : 5]* 





4.66 


0.35 


[0 : 5]** 





3.22 


0.67 


[15 : 20]* 





2.99 


0.78 


[15 : 20]** 





2.36 


0.81 



to the detailed spherical approximation are given, com- 



100 


measured axial bow shock extention 
- 3 Myr fit 
17-20 Myr fit 


= / 


10 


-''"*" / 

...-"+"" 






1 


10 




t [Myr] 





Figure 16. Bow shock extention on the axis versus time. The 
two fits are: 10.45 t ' 49 (0-3 Myr), and 3.49i x 15 (17-20 Myr). 

putcd by application of l|15|) . According to that, the expo- 
nent for the first five million years should be 0.67. Using 
the pure power law, an exponent of 0.35 is achieved in the 
simulation data. Allowing for the radial offset gives a best 
fit exponent of 0.76. Since the exponent of 0.35 is much 
below any expectation, it follows that the initial condition 
is still important in that phase, and the fit with offset is 
more appropriate. The concurrence of the curves increases 
with time and for the last five million years, the exponent 
for the power law fit of the simulation data (0.78) differs 
from that of the spherical approximation by only 0.03. 
From the increasing aspect ratio, an exponent lower than 
the one of the spherical approximation should be expected. 
The simulation shows that the effect is small. 

4.2.5. Axial bow shock propagation 

Accordi ng to self-similar hydrodynamic jet models (com- 
pare e.g. I Alexanderll2002l and references therein), the ve- 
locity of the jet head should be proportional to t'bow.ax = 
£(«+ 2 )/( K +5) ; w here n is the exponent of a local power law 
for the external density (p oc r K ). For a jet head with 
constant area, the velocity should be given by the one- 
dimensional estimate: «bow,ax = v^bcam oc £ 2 /( K + 2 ). The 
last proportionality holds for n > — 2, only. 

The bow shock propagation in axial direction is shown 
in Fig 1161 Its evolution can be represented by two power 
laws that are much closer to the expectation from the 
self-similar models (compare Table OJ than to the one- 
dimensional estimate. For the time range up to three mil- 
lion years, the bow shock is still spherical. The exponent 
should therfore be exactly —0.4. But the initial conditions 
set up the jet with a finite amount of energy from the be- 
ginning which would cause an exponent of —0.6, analogous 
to the supernova case. The measured exponent of —0.51 
shows that the jet is just in transition between these two 
phases. Towards the end of the simulation the jet head 
velocity grows slightly faster than self-similar. The good 
agreement with the self-similar models is probably due 
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Table 3. Comparison between the fitted exponent e for 
the axial bow shock velocity (vbow,ax oc t e ) and the ex- 
ponent expected in the self-similar and one-dimensional 
(ID) estimate, f denotes denotes velocity growth faster 
than any power law. 



time range [Myr] 


measured 


self-similar 


ID 


[0:3] 


-0.51 


-0.4 





[17:20] 


0.15 


0.09 


T 



to the unstable beam, at least in the later phases, where 
the momentum is distributed over a large and increasing 
area. Jet models with intact beams (at higher jet density) 
have been shown to have head advance speeds more con- 
sistent with the on e-dimensional estimate th an with the 
self-similar models l|Carvalho fc O'DealE002a|l . 



4.2.6. Aspect ratio 

Defining the aspect ratio to be the length divided by the 
width of the bow shock, it grows proportionally to t 34 to- 
wards the end of the simulation. This is clearly a non-self- 
similar feature, and may be partially due to the squeezing 
of the cocoon by t he gravity of the swept up gas, as dis- 
cussed recently bv I Alexander! 1 20021) . 



4.2.7. Evolution of the contact discontinuity. 

The temporal evolution of the contact discontinuity is 
shown in Fig. El Both the maximum and the average ra- 
dius arc shown. The position of the contact discontinuity 
was determined by a passive tracer variable that is ad- 
vected with the flow. Due to mixing, the cocoon plasma 
does not keep its initial value of zero. The contact sur- 
face is therfore defined to be at the largest radius with 
the tracer having less than 10% of the ambient medium 
value. The average and maximum values of the contact 
discontinuity evolve in a dissimilar way, which is due to 
the geometry changes discussed above. 

Other than at the earliest times, the contact discon- 
tinuity is always decelerating. For some time (up to t = 
1 Myr) this is sufficient to stabilise against the Raylcigh- 
Taylor instability. To be Raylcigh- Taylor stable, the de- 
celeration of the contact discontinuity must overcome the 
gravity of the assumed dark matter halo, which is given 
by 



5D 



3/?fcT r/a 
/jwiHi 1 + r 2 /a' 



= 1.9410" 



via , o , 

-^cm/s 2 (16) 



l + r 2 /a 2 



The central region is represented by the maximum value, 
which shows the strongest deceleration at the beginning. 
We determine a characteristic deceleration from a fit to 
the maximum value of the contact discontinuity up to 
1 Myr. Note that this is also an upper limit for the 
whole simulation. The result is compared to gravity in 



Fig. EJ This shows that the contact surface should be 
Rayleigh- Taylor unstable for t > 2 Myr, which is con- 
sistent with the density plots. The linear growth time 
(~ ^/l/g = 2.5 Myr ^/(Z/25 pc)(10" 6 cms" 2 /g), where I 
is the wavelength of the instability and g the total accel- 
eration) is typically below the simulation time for wave- 
lengths greater than the resolution limit. 

The details of the acceleration of the contact disconti- 
nuity are more complex than in the global discussion of the 
previous paragraph. Figure IT51 shows that it is constantly 
shaken with values of the order 10 -5 cm s~ 2 . For stabil- 
ity, the global deceleration has to exceed this value, which 
is the case for a fraction of the first Myr only. Another 
complication is that the contact surface is not smooth but 
Kclvin-Hclmholtz fingers always penetrate through. The 
conclusion is that, at least in the central region, swept up 
gas is constantly entrained into the cocoon. 

The entrainment rate is shown in Fig.^5] It is slightly 
rising (oc < 32 ). The entrained mass is compared to the 
mass of the gas that filled the cocoon volume before the 
jet activity. This mass fraction is linearly rising, and would 
reach unity after 1.6 Gyr, if it would continue in the same 
way. 

The width of the cocoon relative to that of the bow 
shock radius (A) is also given in Fig. El For self-similar 
evolution this val ue should be constant, and detailed 
self-s imilar models l|Heinz et alll998tlKaiser fc Alexander! 
Il999l) place it in the range of 0.8 to 0.9. Such a high 
value is reached at no time. The relative width decreases 
with time to 0.4 (0.26) in the case of the maximum (aver- 
ag e) value. Thi s resul t is in agreement with simulations 
by IZanni et all l)2003|) . For the late phase, the average 
value for A levels off. This may indicate a phase of nearly 
self-similar behaviour, only distu rbed by the grav ity of 
the swept up gas, as discussed bv lAlexande^l (|2002|) . The 
slowly growing aspect ratio also supports this view. The 
small value of A may be d ue to the weak bow shock, as 
discussed in more detail bv lZanni et alJ 1 20031) . 



5. Comparison to observations 

5.1. Cygnus A 

This section refers to published data o n Cygnus A 
(|Ca,ri11i Rr. Ba,rthellll99ri fSmith et al.ll2002l) . For illustra- 
tion, an overlay of the smoothed X-ray and radio data is 
shown in Fig. [3U] below the synthetic X-ray image of the 
2.5D simulation. 

The X-ray-radio comparison shows the jct-IGM inter- 
action, impressively. Coincident with the radio cocoon, 
there is a deficit of X-ray emission (disregarding for the 
moment the bright filaments in the center and the emis- 
sion from t he beam and t h e core ) . As discussed above, an d 
detailed in IClarke et all l)l997|) and IZanni et all l)2003|) . 
the parameter £, i.e. the thickness of the shocked ambi- 
ent gas region over the bow shock radius, should exceed 
0.38. Although the 2.5D simulation presented here differs 
from Cygnus A because of the instable beam, it should 
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Figure 17. Position of the contact discontinuity over time. The contact discontinuity is determined by a passive tracer. The 
maximum and the average value is shown. The left figure shows absolute values, the right one shows the value relative to the 
maximum sideways bow shock radius. The fit in the left figure is a power law with exponent 0.28. 
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Figure 18. Left: Comparison of the deceleration (determined by a global fit for the time span up to one Myr) of the contact 
surface against the local gravity (at the average position of the contact surface) of the dark matter halo. Right: Detailed 
acceleration averaged over 0.2 Myr (because of the spatial resolution) for both, the maximum and the average position of the 
contact surface. 



be appropriate for the overall morphology. In the simula- 
tion, the s hocked ambie n t gas h as elliptically shaped X-ray 
isophotes. ISmith et alJ l|2002l) report elliptically shaped 
isophotes from the cocoon boundary out to 66/ (h/0.7) kpc 
(i.e. the outer boundary of their annulus four). Further 
out, the spherical fit is the better one. Identifying the bow 
shock's radius with that position, together with a cocoon 
width of w 32/(/i/0.7) kpc, results in £ = 0.76, in agree- 
ment with the discussed limit. Taking the average value 
for the position of the contact discontinuity, the 2.5D sim- 
ulation produces £ = 0.74 at t — 20 Myr, showing that 
such a high value is indeed possible for sources of that 



size. 



Using the knowledge of the sideways bow shock posi- 
tion, its aspect ratio (length to width) may be derived, 
which turns out to be 1 .1 (using the hot s pot s eparation 
of 147/(ft/0.7) kpc from lCarilli fc Barthell ljl996|) . it is 1.1 
for the easte rn and 1 . 2 for the western jet). From the low 
aspect ratio. iKrausd (|2003l) has concluded a density ratio 



of rj < 10~ 3 . However, also the r\ = 10~ 4 2.5D simula- 
tion presented here reaches an aspect ratio of 1.8 at the 
end of the sim ulation. Since aspect ratios can only grow 
(iKrauseiEiol . and the simulated jets are smaller than 
the observed jet in Cygnus A, it follows that the real jet 
still encounters more mass than the simulated ones. This 
is certainly due to a shallower densit y profile in the at- 
mospherc of Cygnus A of /? = 0.51 (|Smith et alJl2002|) 
compared to (3 — 0.75 employed for the simulations (the 
latter was the old value from the ROSAT data). At the 
largest extent of the bow shock in the 2.5D simulation 
this would increase the ambient density by a factor of 
three, which already might result in the desired reduc- 
tion of the aspect ratio. A still lower central value for the 
density contrast cannot be precluded. However, a density 
contrast of -q ~ 10~ 4 would also b e supported by the large 
cocoon width (|Rosen et al.lll99^) . A nother constrai nt for 
the density contrast was presented bv lKrausel ( 200ril) . who 
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Figure 19. Entrained gas mass in the cocoon over time. Left: Absolute values, the fit function is: 0.677(f/Myr) 1 ' 32 . Right: 
Percentage with respect to the mass of the initial condition in the same volume, the fit function is: 1.5 + 0.06(t/Myr). 



showed that i)>3x 10~ 6 , because the jet is no longer in 
the spherical bow shock phase. 

It has been shown above that the sideways, pressure- 
driven part of the bow shock satisfies with good accu- 
racy even at an aspect ratio of 1.8. Given the low aspect 
ratio of Cygnus A, this equation should be very accurate 
here, and can be used to derive the jet's power and age. 
In order to do that, one needs information on the gas dis- 
tribution before the jet activity. One requirement is that 
the distribution should join the density distribution far 
from the radio lobes, smoothly. Assuming it was a King 
distribution JSJ, as found for other clusters of galaxies, the 
following equation has to be satisfied: 



Pco = 1-64 (a/10 kpej-^mp/cm 



(17) 



The core radius a (for the density profile before the 
jet activity) cannot have been much greater than w 
20/(/i/0.7) kpc. This follows from the enhanced X-ray 
emission next to the radio cocoon. Such enhancements 
have been shown analyt ically to appear in atmospheres 
steeper than p oc r~ - 1 l)Alexa,nderll2nf)3 and references 
therein). In the 2.5D simulation, the line-of-sight inte- 
grated X-ray emission is strongest next to the radio lobes 
for positions of the contact discontinuity in excess of 60% 
of the core radius. There is no obvious lower limit to the 
core radius. However, it will be shown in the following 
that the exact value of a is not needed for the present 
discussion. 

Approximating the gas density p C Q with a constant, 
one infers the following conditions for the jet's power and 
age: 



100tt 



27 



3r bo 
5«bo 



-Pc,O r bow W bow 



(18) 



(19) 



An interesting feature of these equations is that the source 
age is independent of the external density, which remains 



true for the normalisation of arbitrary density distribu- 
tions. Applied to the simulation data, this method yields 
quite accurate values, as demonstrated in sect. 13.4.31 In 
order to apply it to Cygnus A, the bow shock's sideways 
radius is taken to be 66 kpc (see above), and the veloc- 
ity from the temperature jump meas ured from th e X-ra y 
data at that position. According to ISmith et al.1 l)2002|) , 
the temperature of the preshock gas is at about 7-8 keV 
at large radii, falls to 5.4 keV immediately before and then 
rises significantly to 9.2 keV at the position where the bow 
shock was proposed to be situated above. It is therefore a 
wea k shock. Th e shock conditions for a weak shock yield 
(e.g.[Sh3[l9923): 
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The measurement indicates a temperature jump T\/Tq of 
1.7, and errors limit it to smaller than 4. Also, the bow 
shock velocity should exceed the speed of sound in the 
preshock gas: c s = 1313\/Xb/5.4keV km/s. This results 
in a sideways bow shock velocity of: 



2217 ±Hf km/s 



Inserting this into Ijl8(l and l|19|l and assuming an average 
pre-jet density of p c ,o = 0.05m p /cm 3 yields L = 4.4 ±3 5 5 
xlO 47 erg/s and t = 17=F? 2 Myr. 

The parameters can also be estimated assuming a King 
distribution in the pre-jet era. It follows from Q: 



L 



47r PcO r bow^bow I(r how /a) 3 
9a J(rbow/a) 2 
3a 2 J(r bo w/a) 
^bow^bow I(rbow/a) 

" x 2 (l + x 2 )°- 765 dx 
o 



(20) 
(21) 
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J(z) 



where p e ,o is given by i|17|) . The result is shown graphically 
in Fig. [21 It turns out that all the parameters are only 
weakly dependent on the core radius. The average value 
for each curve is given in the following summary: 



44.3 
6.3 

16 Myr 



xl0 47 erg/s 



L = 7.9 ± 

t = 24=F n 



E = 5.9 ± 44 8 9 10 62 erg 



where E is the total energy released by the jets. As dis- 
cussed above, in order to obtain the power of the beam, 
the flux of internal energy through the bow shock has to be 
subtracted and the gravitational energy increase has to be 
added. A reasonable estimate of the internal energy con- 
tained within the bow shock region before the jet activity 
can be obtained by assuming an isothermal King atmo- 
sphere. Adopting a core radius of 10 kpc, the X-ray data 
far from the center is consistent with a central pressure 
of po = 10~ 9 erg/cm 3 . This leads to a total of 10 60 erg, 
which is negligible here. The gravitational energy increase 
is not so easy to estimate. But for the present discussion it 
is sufficient to say that the 2.5D simulation indicates that 
up to 30% of the jet power is used to lift up gas. Adopting 
a value of 20%, the power through the jet channel is de- 
termined to L > 1.9 x 10 47 crg/s. Neglecting this rather 
uncertain contribution, the lower limit to Cygnus A's jet 
power is 1.6 x 10 47 erg/s. 

This has to be compared to the power 
in the jet channel. For a non-relativistic jet: 



5 x 10 45 (i?j/0.57kpc) 2 (10 4 ?7o)(pc,o/m F 
^/(c/2)) 3 erg/s . 



,-3\ 



The j et's radius was adopted from ICarilli fc Barthell 
{[1996^ usin g the latest value for the Hubble constant, 
h = 0.72 llSnergel et a.l . 1 l200.3h . and the subscript zero 
denotes values in the center. Hence, even when using 
most optimistic numbers, the kinetic jet power falls short 
of the lower limit derived above, by a factor of ten or 
more. This means that the internal energy cannot be 
responsible for the rest, since this would yield a subsonic 
jet. Therefore the energy should be in the (toroidal) 
magnetic field: 



L, 



2tt RfuBT'vj 



5 x 10 44 (i?j/0.57kpc) 2 /?r 2 (u B /u BiHS ) 



Here, T is the bulk Lorentz factor f3 = Wj/c, and Ub HS = 
9x 10 -10 erg/cm 3 is t he magnetic e n ergy d ensity measured 
in the hot spots by I Wilson et alJ l|2000|) . It is again an 
upper limit, since the magnetic field in the jet will be 
lower than in the hot spot, where it is shock compressed. 
It is clear that the derived power can also not be raised 
by the magnetic energy flux, as long as we assume a non- 
relativistic jet. 
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Figure 21. Detailed results for age (top), and power (bottom) 
for Cygnus A's jets, assuming h = 0.7 and the sideways bow 
shock velocity indicated on the figures. 



Therefore, the jets in Cygnus A have to be relativistic. 
The kinetic power for a relativistic jet may be written as 
follows. 



L 



kin.R — 



27nRjWe(l - (r/i)- x )/3c 3 
8.8 x 10 46 (i?j/0.55kpc) 2 (10V, ) 
(p e ,o/m p cm- 3 )(l - (r/i)- 1 )/? erg/s 
PjhT 2 / Pc , 



where h = 1 + je/pc 2 is the specific enthalpy, t^r de- 
notes the relativistic generalisation of the density contrast. 
This number should be similar to the valu e for the non- 
relativistic r\ (compare iRosen et al.lll999|) . The central 
cluster density p .o has a weak dependency on the core 
radius which is subsumed in a variation of 7/r.o, which is 
regarded to have a value below 10~ 4 . This can be com- 
bined with constraints on the relativistic Mach number 
M = r/3/r s /? s , wher e the index s den otes the values for 
the sound speed fe.g. lMarti et al.l|l997tk The sound speed 
is given by c 2 = jp/hp = (7 — 1)(1 — l/h)c 2 , which yields 
for the specific enthalpy: 

M 2 + T 2 /3 2 ™ - "> 

h = TT7T- ^7^77, 9 



M 2 + gT 2 f3 2 ' 



7-1 



(22) 
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ate one, showing huge radio lobes like FR II sources, but 
no edge brightening. It seems to be in a similar state to the 
2.5D simulation towards the end, presented above. The jet 
is unstable, probably because it cannot reach the required 
Mach number in order to stabilise it up to the tip of the 
cocoon. The whole cocoon seems to be turbulent, and en- 
training external gas across the contact discontinuity. This 
might be the reason for the lack of an X-ray deficit as ob- 
served in Cygnus A. However, the ROSAT data may lack 
the necessary resolution. What it does show is that the 
X-ray gas is elongated in the direction of the radio source. 
This can be interpreted as a wide region of shocked ambi- 
ent gas, as in Cygnus A. 



Figure 22. Power through the jet channel of Cygnus A over 
the Lorentz factor compared to the lower limit derived in the 
text. Curves for differing Mach numbers are almost identical. 



Since Cygnus A shows a stable and laminar jet, the Mach 
number has to exce ed unity. Most likely the Mach number 
is greater than that. ICarilli fc Barthell jl99rfl infer a Mach 
number of M w 8 (from the oblique shocks) and M < 13 
(from the opening angle). The total power derived in this 
way is plotted against the Lorentz factor in Fig. [221 f° r 
different values of the central density contrast and Mach 
number. The lower power limit derived above results in a 
Lorentz factor around 10, the central value gives 37 to 39. 
An important result is that the jet's Lorentz factor should 
exceed ten for reasonable assumptions. According to (|22|l 
M should then exceed 14, in order to get h > 1, which 
might point to a significant Alfven speed. 

If one would regard a Lorentz factor above twenty as 
unreasonable, the jet power could be further constrained 
to less than 3 x 10 47 erg/s. 

The derived jet power exceeds the total emit- 
ted radio power by at least a factor of w 100 
I Carilli fc Barthell Il99fil). and might be compatible with 
esti mates bvlKaiser fcAlexandeilll 19991) (2 x 10 46 erg/s) 
and lZanni et all l|2003|) (3 x 10 46 erg/s). who do not state 
their errors. 

The simulations have shown that the Kelvin- 
Hclmholtz instability produces large fingers in the central 
regions of the cocoon. These may be identified with the 
belts of cooler (4 keV) gas inside the cocoon of Cygnus A. 
Additionally, the contact discontin uity might be Rayleigh- 
Taylor unstable (| Alexander! l2002f) . The 2.5D-simulation 
shows that the deceleration is stronger at the beginning, 
and falls below the critical value for most of the simulation 
time. Also, local shaking dominates over the global decel- 
eration. Therefore the contact discontinuity is expected to 
be Rayleigh- Taylor unstable. 



5.2. Hercules A 

Radio and ROSAT data o n Hercules A are given by 
iGizani fc Leahvl 1)20031120041) . This source is an intermedi- 



al Abell 2052 

X-ray, radio, emission line, a nd other data on this centra l 
cluster galaxy can be found in lBlanton et alJ ||2001LEoQ3]) . 
It is an old radio galaxy, the radio emission is diffuse, and 
no sign of hot spots or beams has been found. The diffuse 
radio emission is d evoid of X-rays, and is sur rounded by 
bright X-ray shells. iBlanton et all l|200lll2003|) argue that 
the shells consist of two oppositely situated bubbles, as 
found in other X-ray clusters. However, the study of the 
dependency of the X-ray emission on the viewing angle 
suggests that the situation may be similar to the ten de- 
gree plot of Fig. El A low density jet may have blown an 
approximately spherical bubble. After reaching the criti- 
cal radius (see above), a jet with a narrow beam, i.e. high 
thrust bored a cigar-like extension into the shell. Now, 
we view the cigar-like extension from such an angle that 
it partly overlaps with the big shell, producing enhanced 
emission at those parts of the X-ray rings. 

An interesting point in this interpretation is that when 
modeling the X-ray enhanced regions by two elliptical, 
overlapping rings, these are elongated in the same direc- 
tion, and deprojecting to a circular ring results in an in- 
clination of roughly 37° for both. This inclination is com- 
patible with the estim ate of « 45° from the radio data by 
IBlanton et all l|2003|) . 

A question concerning this interpretation arises from 
the temperature map of the X-ray shells. The smaller shell 
seems to be colder than the larger one. While the uniform 
temperature of each of the elliptical shells is consistent 
with their proposed nature, the cigar part, correspond- 
ing to the smaller elliptical ring, is always hotter in the 
simulations than the inner bow shock part. This problem 
might be alleviated by the fact that the jet in this source 
is presently not active. Therefore, the radio plasma that 
was present in the cigar may have left it in the backflow, 
leaving behind an undcrprcssured region. The gas in the 
shell would then expand somewhat into the cigar cavity, 
thereby lowering its pressure. For that reason, the expan- 
sion velocity of the leading bow shock would slow down. 
Assuming a cylindrical shell for the cigar part, the adia- 
batic expansion (T oc V~ 8 / 3 ) would require a change of 
the inner radius of the shell from e.g. 80% to 70% of the 
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outer one, in order to explain a temperature change by a 
factor of two. 

Following the smaller X-ray shell that has been inter- 
preted as cigar part above , there are emission line fila- 
ments i Blanton et al-lEoOll) . These might be due to weak 
radiative shock waves inside of the shell that might have 
been excited by the re-expansion described above. 

5.4. Radio galaxies at higher redshift 

Radio galaxies with redshifts greater than z w 0.6 show 
associated large sc ale emission line r egions that are aligned 
with the radio jets llMcCarthvll993l) . The aligned emission 
line regions are often located where one would expect the 
radio cocoon, an d are stronger towards the center (e.g. 
iBest et al"1ll996j) . It has been shown in this paper that 
the central cocoon regions can entrain several percent of 
the gas mass that was situated there before the jet event. 
The X-ray belts in Cygnus A have been interpreted above 
as evidence for this entrainment process. These belts are 
the coldest X-ray gas found in the Cygnus A system. It 
is therefore possible that such entrained gas reaches ther- 
mal instability, if the cluster's gas density is higher, or the 
source has some more time to cool. The cooling time for 
these belts is some 100 Myr. One can therefore imagine a 
situation where the entrained gas cools down to emission 
line temperatures whereas the other shocked ambient gas 
does not cool. As suspected above for Abell 2052, the cool- 
ing could also proceed within X-ray filaments that con- 
tain weak, radiative shocks. This would explain why the 
line ratios for some sources (preferenti ally the smaller) ar e 
consistent with shock excitation (e.g. Ilnskip et al.ll2002h . 
without relying on the compression of pre-existing clouds 
by the bow shock. Wether this picture can cope with other 
observational data has yet to be explored. 

6. Summary and conclusions 

Bipolar jet simulations have been presented in 
2.5D and 3D. For the 3D simulation, a cylindrical 
grid was employed, and the boundary conditions were 
given. It has been shown that the influence of the axial 
boundary remains acceptably small. Because of small 
disturbances in the ambient medium, the backflows are 
located at different distances from the jet axis, and 
permeate each other in the center. There, a turbulent 
mixing region emerges that entrains shocked ambient gas 
across the contact discontinuity via Kelvin-Hclmholtz 
and Raylcigh- Taylor instabilities. After a certain time, 
both simulations show two bow shock phases. The shape 
of the outer one is well fitted by a parabola. The inner 
part has an elliptical shape and stays axi-symmetric, even 
in the 3D-simulation. It follows the blastwave equation 
of motion with good accuracy (a few percent difference 
in the exponent of a local power law fit), even when 
the aspect ratio is high. The aspect ratio keeps growing 
for the whole simulation time, i.e. a self-similar regime 
is not reached. The viewing angle-dependent emission 



maps show that the different parts of the bow shock may 
form circular and elliptical rings, and produce elliptical 
isophotes when viewed at high inclination. Varying the 
inclination can also produce X-ray deficits. At late times 
of the 2.5D simulation, the beam is unstable, barely 
reaching the tip of the bow shock. This could be expected 
from stability analysis, because a non-relativistic very 
light jet cannot keep a high Mach number, which is 
necessary for stability. Regarding the propagation of the 
tip of the bow shock, an armlength asymmetry of a few 
percent was measured in both simulations. For late times, 
the jets are faster on the side with the on average higher 
density. 

The X-ray structure in Abell 2052 was shown to be ex- 
plainable by two elliptical rings, one from the outer bow 
shock part of a former jet and one from the inner one. 
From the ellipticity, an inclination of 37° is derived con- 
sistent with other estimates in the literature. Based on the 
simulation results, the low inclination is one reason for the 
pronounced X-ray deficit inside of the rings. 

The X-ray gas in Hercules A and Cygnus A shows el- 
liptical isophotes elongated in the direction of the radio 
jets. The ellipticity decreases at a certain distance from 
the center. The location where the ellipticity drops was 
identified as being the bow shock position. For Cygnus A, 
the quality of the data is sufficient to determine the side- 
ways bow shock position to 66 kpc. At the same position, a 
temperature jump can be infercd from the X-ray data, im- 
plying a sideways bow shock velocity between 1300 km/s 
and 4200 km/s, i.e. a Mach number between one and three. 
The width of the radio cocoon is only a quarter of the bow 
shock width, which is consistent with the observation of an 
X-ray deficit, and shown to be possible in the 2.5 D simu- 
lation presented here. However, in self-similar models, this 
fraction is usually above 80%, which is a shortcoming of 
these models. Applying the simulation results and using 
the width and velocity of the sideways bow shock, a jet 
power of > 10 47 erg/s and an age of sa 24 Myr is derived. 
This result was found to be inconsistent with the jet being 
non-relativistic, and a lower limit for the Lorcntz factor of 
ten was infered. 

Explaining the belts of low temperature X-ray gas 
within the radio cocoon of Cygnus A by entrained gas over 
the contact discontinuity, it was suggested that such gas 
could cool further to form emission line regions, if the den- 
sity of the ambient gas would be somewhat higher. This 
scenario could be realised at higher redshift and explain 
the origin of the gas producing the alignment effect. 
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Appendix A: Late phase of spherical blastwaves 

For which radii is it possible to use the blastwave ap- 
proximation, i.e. the solutions of This equation was 
derived under the assumption of vanishing external pres- 
sure, which is equal to a strongly supersonic shock. In the 
simulations we get weak bow shock Mach numbers very 
soon. Then, the external pressure (p e xt) and possibly also 
the gravitation (remember that the atmosphere is held in 
equilibrium of pressure and gravity) has to be taken into 
account. This will be discussed in the following for the 
case of isothermal power law atmospheres, applicable to 
the presented work. 

The force ballance equation for a shell, driven by en- 
ergy input E{t) = Ct d into an environment with gas mass 
profile M(r) = 4nr 2 p(r)dr, where the density distri- 
bution is given by p(r) = po(r/ro) K , can be written as 
follows: 



d 

— (Mv) = S(p int -Pext) 

at 



GMbmM 



(Al) 



Here, S — Airr 2 , G is the gravitational constant, A^dm 
denotes the gravitating Mass profile of the dark mat- 
ter halo, and the internal pressure is given by p lnt = 
2(E(t) - Mv 2 /2)/W, where V = 4irr 3 /3. Using hydro- 
static equilibrium, this equation can be nondimension- 
alised and rearranged in the following way: 
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The scales ro and to are given by: 



r = ^J%k + 5)/5c s i 
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For typical values of galaxy clusters, say po = 
10~ 26 g/cm 3 , constant sound speed of c s = 1000 km/s, 
and C = 10 47 erg/s for d — 1 the typical scale is Mpc, 
reducing by a factor of a few when the jet turns off and 
the energy stays constant at some 10 62 erg (d = 0). The 
force balance equation was numerically integrated, an ex- 
ample is shown in Fig. IA.1I The general result is that 
blast waves with finite energy fall back after reaching the 
critical radius, the formal solution being oscillatory. The 
ones with constant energy injection break and change the 
slope of their power law at that radius. Up to that radius, 
the propagation follows the same power law that can be 
derived using the strong bow shock hypothesis. The con- 
clusion is that for the simulated jets, and the typical real 
sources, the critical radius has not yet been reached, and 
the blastwave equation of motion can be expected to be 
a good approximation as long as non-spherical effects re- 
main small. 
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Figure A. 1. Numerical solutions of the force balance equa- 
tion for late spherical blast waves. The blast wave with con- 
stant energy injection changes the power law slope at the crit- 
ical radius. The one with finite energy falls back. 



